We’ll use the following data and R packages today:
https://github.com/lekoenig/Rspatial-UNH.git
library(sf) library(dplyr) library(ggplot2) library(ggspatial) library(mapview) library(rnoaa)
11/15/2021
We’ll use the following data and R packages today:
https://github.com/lekoenig/Rspatial-UNH.git
library(sf) library(dplyr) library(ggplot2) library(ggspatial) library(mapview) library(rnoaa)
1. geospatial analysis and operations
1. geospatial analysis and operations
2. visualizing spatial data
Example field site: Swains Lake in Barrington, NH
Example field site: Swains Lake in Barrington, NH
Example field site: Swains Lake in Barrington, NH
Load {sf} and read in a subset of the NH National Hydrography Dataset (available from www.granit.unh.edu):
#install.packages("sf")
library(sf)
local_lakes <- st_read("./data/NHDWaterbodyDataset_NHSeacoast.shp")
Load {sf} and read in a subset of the NH National Hydrography Dataset (available from www.granit.unh.edu):
#install.packages("sf")
library(sf)
local_lakes <- st_read("./data/NHDWaterbodyDataset_NHSeacoast.shp")
class(local_lakes)
Load {sf} and read in a subset of the NH National Hydrography Dataset (available from www.granit.unh.edu):
class(local_lakes)
## [1] "sf" "data.frame"
Load {sf} and read in a subset of the NH National Hydrography Dataset (available from www.granit.unh.edu):
head(local_lakes)
Load {sf} and read in a subset of the NH National Hydrography Dataset (available from www.granit.unh.edu):
head(local_lakes,3)
## Simple feature collection with 3 features and 13 fields ## Geometry type: POLYGON ## Dimension: XY ## Bounding box: xmin: -71.28925 ymin: 42.91414 xmax: -70.7096 ymax: 43.44024 ## Geodetic CRS: NAD83 ## Prmnn_I FDate Resoltn GNIS_ID GNIS_Nm AreSqKm Elevatn ## 1 141035025 2005-10-04 2 <NA> <NA> 0.00100000 NA ## 2 141036768 2019-03-19 2 <NA> <NA> 0.00028992 NA ## 3 141032014 2019-03-19 2 <NA> <NA> 0.00038662 NA ## ReachCd FType FCode VsbltyF Shp_Lng Shap_Ar ## 1 01060003011689 390 39004 24000 0.0010302955 7.489623e-08 ## 2 <NA> 466 46600 0 0.0007273591 3.203682e-08 ## 3 01060003008954 390 39004 50000 0.0007978068 4.298463e-08 ## geometry ## 1 POLYGON ((-70.86617 42.9143... ## 2 POLYGON ((-71.28921 43.0608... ## 3 POLYGON ((-70.70971 43.4399...
Interactive mapping tools are great for sharing and communicating spatial data, but also when you just want to look at the data.
One of these packages is
{mapview}. Let’s take a look at our local_lakes dataset:
#install.packages("mapview")
library(mapview)
mapview(local_lakes)
{sf} integrates well with ‘tidy’ data workflows. What do I mean by that?
{sf} integrates well with ‘tidy’ data workflows. What do I mean by that?
Let’s subset our local_lakes dataset to only include my field site, Swains Lake:
library(dplyr) swains_lake <- local_lakes %>% filter(GNIS_Nm == "Swains Lake")
{sf} integrates well with ‘tidy’ data workflows. What do I mean by that?
Let’s subset our local_lakes dataset to only include my field site, Swains Lake:
library(dplyr) swains_lake <- local_lakes %>% filter(GNIS_Nm == "Swains Lake") mapview(swains_lake)
The st_read() function automatically stores information about our data, including:
The st_read() function automatically stores information about our data, including:
st_geometry_type(swains_lake)
## [1] POLYGON ## 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE
The st_read() function automatically stores information about our data, including:
st_crs(swains_lake)
## Coordinate Reference System: ## User input: NAD83 ## wkt: ## GEOGCRS["NAD83", ## DATUM["North American Datum 1983", ## ELLIPSOID["GRS 1980",6378137,298.257222101, ## LENGTHUNIT["metre",1]]], ## PRIMEM["Greenwich",0, ## ANGLEUNIT["degree",0.0174532925199433]], ## CS[ellipsoidal,2], ## AXIS["latitude",north, ## ORDER[1], ## ANGLEUNIT["degree",0.0174532925199433]], ## AXIS["longitude",east, ## ORDER[2], ## ANGLEUNIT["degree",0.0174532925199433]], ## ID["EPSG",4269]]
Proper representation of coordinate reference systems (CRS) is very important for all spatial work!
CRS are often denoted using an “EPSG” code (see bottom of st_crs output). Common CRS include:
See other resources for more information, including https://datacarpentry.org/organization-geospatial/03-crs/
st_crs is your friend.
st_crs(swains_lake)
The st_read() function automatically stores information about our data, including:
st_bbox(swains_lake)
## xmin ymin xmax ymax ## -71.05927 43.18261 -71.02456 43.19972
The st_read() function automatically stores information about our data, including:
swains_lake
## Simple feature collection with 1 feature and 13 fields ## Geometry type: POLYGON ## Dimension: XY ## Bounding box: xmin: -71.05927 ymin: 43.18261 xmax: -71.02456 ymax: 43.19972 ## Geodetic CRS: NAD83 ## Prmnn_I FDate Resoltn GNIS_ID GNIS_Nm AreSqKm Elevatn ## 1 141033653 2015-01-14 2 00870298 Swains Lake 1.379346 NA ## ReachCd FType FCode VsbltyF Shp_Lng Shap_Ar ## 1 01060003002342 436 43600 1000000 0.1605085 0.0001527365 ## geometry ## 1 POLYGON ((-71.05117 43.1887...
Say we want to find the closest weather station to Swains Lake.
Say we want to find the closest weather station to Swains Lake.
We can use the {rnoaa} package to identify NOAA climate stations, and then use spatial joins to select the nearest site and download some data.
library(rnoaa)
# Takes a minute to run:
stations <- rnoaa::ghcnd_stations() %>%
filter(state=="NH",element=="PRCP",last_year>2020) %>%
st_as_sf(.,coords=c("longitude","latitude"),crs=4326) %>%
arrange(id)
stations <- read_sf("./data/NH_GHCND_stations.shp")
Say we want to find the closest weather station to Swains Lake.
We can use the {rnoaa} package to identify NOAA climate stations, and then use spatial joins to select the nearest site and download some data.
head(stations,3)
## Simple feature collection with 3 features and 9 fields ## Geometry type: POINT ## Dimension: XY ## Bounding box: xmin: -71.55 ymin: 43.4356 xmax: -71.3231 ymax: 43.5541 ## Geodetic CRS: WGS 84 ## # A tibble: 3 × 10 ## id elevation state name gsn_flag wmo_id element first_year last_year ## <chr> <dbl> <chr> <chr> <chr> <chr> <chr> <int> <int> ## 1 US1NHBK0001 207. NH TILT… <NA> <NA> PRCP 2009 2021 ## 2 US1NHBK0002 214 NH BELM… <NA> <NA> PRCP 2009 2021 ## 3 US1NHBK0007 172. NH LACO… <NA> <NA> PRCP 2009 2021 ## # … with 1 more variable: geometry <POINT [°]>
Calculate the distances between each weather station and Swains Lake.
st_distance(stations,swains_lake)
Calculate the distances between each weather station and Swains Lake.
We need to transform the stations and swains_lake data so that they share the same CRS.
# Transform data to CONUS Albers Equal Area Conic projection: swains_lake_proj <- st_transform(swains_lake,5070) stations_proj <- st_transform(stations, 5070)
Calculate the distances between each weather station and Swains Lake.
stations_proj2 <- stations_proj %>% mutate(dist = st_distance(.,swains_lake_proj,by_element = TRUE)) stations_proj2 %>% arrange(dist) %>% select(id,name,dist) %>% head(3)
Spatial join: Find the closest weather station to Swains Lake.
closest_station <- st_join(swains_lake_proj,stations_proj,join=st_nearest_feature)
Spatial join: Find the closest weather station to Swains Lake.
closest_station <- st_join(swains_lake_proj,stations_proj,join=st_nearest_feature) head(closest_station)
## Simple feature collection with 1 feature and 22 fields ## Geometry type: POLYGON ## Dimension: XY ## Bounding box: xmin: 1993730 ymin: 2506799 xmax: 1996588 ymax: 2508389 ## Projected CRS: NAD83 / Conus Albers ## Prmnn_I FDate Resoltn GNIS_ID GNIS_Nm AreSqKm Elevatn ## 1 141033653 2015-01-14 2 00870298 Swains Lake 1.379346 NA ## ReachCd FType FCode VsbltyF Shp_Lng Shap_Ar id ## 1 01060003002342 436 43600 1000000 0.1605085 0.0001527365 US1NHST0040 ## elevation state name gsn_flag wmo_id element first_year last_year ## 1 80.8 NH BARRINGTON 3.2 E <NA> <NA> PRCP 2016 2021 ## geometry ## 1 POLYGON ((1994469 2507096, ...
Now we can download and plot our precipitation data:
prcp_data <- rnoaa::ghcnd_search(closest_station$id, var = "PRCP") # Plot data: prcp_data$prcp %>% filter(date>"2021-08-01" & date<"2021-10-01") %>% mutate(precip_mm = prcp/10) %>% ggplot() + geom_line(aes(x=date,y=precip_mm)) + labs(x = "Date (2021)", y = "Precipitation (mm)") + theme_classic()
Useful packages for mapping in R:
We can still use {ggplot2} for plotting spatial data.
ggplot() + geom_sf(data=swains_lake) + coord_sf()
We can still use {ggplot2} for plotting spatial data.
A lot of the formatting should already be familiar if you’re a {ggplot2} user:
ggplot() +
geom_sf(data=swains_lake,fill="cyan3",color="cyan4") +
coord_sf(ylim = c(43.18, 43.2)) +
ggtitle("Swains Lake, NH") +
theme_classic()
Let’s load one more spatial dataset to include in our map - a shapefile of the watershed boundary, that is, the area of land that drains into Swains Lake.
swains_watershed <- read_sf("./data/Swains_watershed_boundary.shp")
st_crs(swains_watershed)
## Coordinate Reference System: ## User input: WGS 84 ## wkt: ## GEOGCRS["WGS 84", ## DATUM["World Geodetic System 1984", ## ELLIPSOID["WGS 84",6378137,298.257223563, ## LENGTHUNIT["metre",1]]], ## PRIMEM["Greenwich",0, ## ANGLEUNIT["degree",0.0174532925199433]], ## CS[ellipsoidal,2], ## AXIS["latitude",north, ## ORDER[1], ## ANGLEUNIT["degree",0.0174532925199433]], ## AXIS["longitude",east, ## ORDER[2], ## ANGLEUNIT["degree",0.0174532925199433]], ## ID["EPSG",4326]]
Let’s load one more spatial dataset to include in our map - a shapefile of the watershed boundary, that is, the area of land that drains into Swains Lake.
ggplot() +
geom_sf(data=swains_lake,fill="cyan3",color="cyan4") +
geom_sf(data=swains_watershed, fill=NA, color="gray20", lwd=0.6) +
coord_sf() +
ggtitle("Swains Lake, NH") +
theme_classic()
Add scales, orientation arrows, etc. using {ggspatial}.
You can probably think of other ways you want to edit your map for a report or publication, but we’re definitely getting closer!
library(ggspatial)
swains_map <- ggplot() +
annotation_map_tile(zoom = 15) +
geom_sf(data=swains_lake,fill="cyan3",color="cyan4") +
geom_sf(data=swains_watershed, fill=NA, color="gray20", lwd=0.6) +
ggtitle("Swains Lake, NH") +
theme_classic() +
labs(x="Longitude",y="Latitude") +
coord_sf() +
annotation_north_arrow(location="bl",
height = unit(.75, "cm"),
width = unit(.5, "cm")) +
annotation_scale(location="bl",
pad_x = unit(1, "cm"),
text_cex = 0.5,
height = unit(0.2, "cm"))
Add scales, orientation arrows, etc. using {ggspatial}.
You can probably think of other ways you want to edit your map for a report or publication, but we’re definitely getting closer!